Dynamical Effective Medium Theory for Quantum Spins and Multipoles 
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A dynamical effective medium theory is presented for quantum spins and higher multipoles such 
as quadrupole moments. The theory is a generalization of the spherical model approximation for 
the Ising model, and is accurate up to 0(l/z n ) where z n is the number of interacting neighbors. 
The polarization function is optimized under the condition that it be diagonal in site indices. With 
use of auxiliary fields and path integrals, the theory is flexibly applied to quantum spins and higher 
multipoles with many interacting neighbors. A Kondo-type screening of each spin is proposed for 
systems with extreme quantum fluctuations but without conduction electrons. 



I. INTRODUCTION 



Strongly interacting localized electrons have been attracting renewed interest. A key feature characterizing the new 
development is the coupling of spin and orbital degrees of freedom. For example, Ce^^La-rBg has both magnetic 
and electric quadrupole orders, and the phase diagram exhibits intriguing systematics as x varies. The specific heat 
above the quadrupolar ordering temperature shows a large contribution of fluctuation, which is suppressed by applied 
magnetic field. || Similar but non-identical behavior has been found in TmTe. 0] We note that Cei-^La^Be shows the 
Kondo effect in the resistivity, while TmTc is insulating. In 3d systems, the orbital degrees of freedom appears as the 
;> static and dynamic Jahn- Teller effects. [|J An exemplary system is Lai-^Sr^MnOa, which shows the drastic change of 
I/") ' magnetic anisotropy with increasing x. The transport property is characterized by the colossal magneto-resistance. 
I""""- All these phenomena require simultaneous account of orbital and spin degrees of freedom. 

These developments motivate construction of a new quantum theory which can deal with fluctuation effects of not 
only spin but higher-order multipoles. Concerning previous efforts toward the dynamical theory, we refer to the work 
of Hubbard || which addresses the high-temperature limit of the Heisenberg model in high dimensions. For the 
low-dimensional Heisenberg model, highly sophisticated theories are available with account of quantum fluctuations. 
H However, these theories are not suitable to three-dimensional systems with a different character of fluctuations. 

In the case of strongly correlated itinerant electrons, microscopic theories have been developed with account of both 
the Kondo effect and the intersite correlation for the Anderson lattice, |l0|-|l3[] and for the Hubbard model. |l3|^15|] 
I ■ These theories use the idea of dynamical effective medium which is justified in the limit of large spatial dimensions. 
Although the approach is highly successful in deriving the density of states of electrons with nontrivial structure 
around the Fermi level, the magnetic property is treated rather crudely. To clarify the need of improvement, let us 
take an example of the half-filled Hubbard model. With large Coulomb repulsion the single-particle spectrum has 
an energy gap. As a result the spectrum of the effective medium to determine the Green function also has a gap. 
Under this condition the infinite-dimensional theory predicts a gapful spectrum also for spin excitations. In reality, 
however, the spin excitation is gapless in many cases even though the single-particle spectrum has a gap. Furthermore 
in the insulating paramagnetic phase the entropy derived by the theory does not vanish at zero temperature. With 
a magnetic ordering, the entropy can vanish by the same mechanism as the band antiferromagnetism. However, in 
the limit of vanishing charge fluctuation, one has complete spin polarization in the ground state. This last feature 
of the theory is not satisfying since quantum spin fluctuations should reduce more or less the magnitude of ordered 
moments. In order to remedy the situation one has to go beyond the lowest-order self-consistent theory. Then there 
appears a difficulty of analyticity in the Hubbard model and other fermion models. Jl6[ 

The purpose of this paper is to present a next-leading order self-consistent theory which is free from the above 
deficiency. The formulation has so far been attained only for systems with localized electrons such as quantum spins. 
The present theory can deal with a possible paramagnetic ground state without residual entropy, provided that the 
spin excitation has a gapless spectrum. The paper is organized as follows: In the next section we take the Ising model 
as the simplest system to apply the formalism. We exploit the variational character of the thermodynamic potential in 
deriving the magnetization and the susceptibility. The same variational property is also used to compute the entropy 
and specific heat within the same scheme. Section 3 describes the extension of the formalism to quantum models. 
We take the Heisenberg model with arbitrary exchange interactions as a representative quantum model. The static 
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approximation is introduced which is much simpler to treat than the fully dynamical counterpart. The approximation 
is justified in the high temperature limit and makes a correspondence to the classical theory. In §4 we discuss a way to 
solve the problem by mapping to the spin-boson system where an impurity spin interacts with bosonic environment. 
The latter represents the dynamical medium provided by surrounding spins. By this mapping, numerical methods such 
as the Monte Carlo technique, numerical renormalization group, or the self-consistent perturbation theory become 
applicable to solve the impurity problem. Section 5 applies the theory to more general systems with quadrupole 
moments or with crystalline-electric-ficld effects. In the final section we discuss the results, in particular a possibility 
of a Kondo-type effect in spin systems without conduction electrons. 



II. SELF-CONSISTENT RENORMALIZATION FOR THE ISING MODEL 



A. Perturbation theory with auxiliary fields 



Although our principal interest is in quantum models, we first take the Ising model to present the new formalism 
in the simplest, and hence the most transparent manner. The spherical model approximation (SMA) goes one step 
beyond the mean-field theory for the Ising model, and has been discussed for many years. Jl7| We rederive the SMA in 
a variational formalism, which makes clear the structure as well as limitation of the approximation. As a byproduct, 
we derive a formula to express the specific heat in terms of the susceptibility. 

The Hamiltonian of the Ising model is given by 



(i) 



where each Oi takes ±1 and hi denotes a magnetic field at site i. For clarity we assume that the interactions Jy are 
all positive, and the inverse of the matrix {Jij} exists. The self-consistency equations, however, are valid for more 
general interactions. We introduce auxiliary fields which obey the Gaussian distribution, and which couple with a 
locally. The basic identity is written as 



° xp ( f x j i3Wi I = ( det i 3 j ) 1/2 n f 



exp 



■?5> 



(2) 



where the the integration path C runs from — _Rexp(i7r/4) to _Rexp(i7r/4) with R going to infinity. This choice of a 
path makes it possible to obtain the convergent integral after the change of integration variables. Namely we distort 
the path either along the whole real axis or the imaginary axis depending on the sign of the eigenvalue of the matrix 
{Jij}- The partition function is written as 



(3) 



where 



(4) 



Thus we see that the Gaussian identity casts the original two-body interaction into the sum of local coupling terms. 

In order to perform perturbation theory we first set the energy scale to which corresponds to the ferromag- 

netic transition temperature in the mean- field approximation (MFA). Thus each is of order l/z n where z n is the 
number of equivalent neighbors. To organize the perturbation series concisely, it is convenient to use the fermion 
representation for Ising spins. Namely we introduce spinless fermions by 



a i = 2flf i -l, 



(5) 



where fi is the annihilation operator of a fermion at site i. Then each term of the perturbation expansion can be 
expressed in terms of Feynman diagrams. The intersite contribution involves the bare cumulant {4'i4 > j)o = TJij- We 
remark that the same diagrams appear when we use the fermion representation from the outset without introducing 
the <f> field. However for quantum models the cf> field makes the theory much more concise. 
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Figure 1(a) shows the zero-th order contribution to the polarization function which corresponds to the self-energy 
of the renormalized exchange interaction. The small parameter l/z n is offset by the site summation in these diagrams 
which are called tadpoles. The tadpoles lead to shift of the effective fermion level, and represent as such the molecular 
field. 
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FIG. 1. Feynman diagrams for the polarization function, (a) Contributions of the zeroth order in l/z„; (b) first order; (c) 
second order. The solid lines represent the Green function of fermions, and the dashed lines do the exchange interaction. 

The set of zero-th order diagrams make a tree- like structure composed of tadpoles. In the leading order there 
are no loops composed of interaction lines. Summation of all the tree diagrams is equivalent to the saddle point 
approximation of eq.(||). Namely after taking the trace over ui in eq.(f|), we find that the saddle point of the field 
4>i satisfies the condition 

- tanh^Oi + a,) = 0. (6) 

3 

It is seen that a,i gives the molecular field at site i, and the magnetization m, = (ui) is given by nii — tanh/3(/i.; + a^). 
Namely we get 

•i, y "] ./, ; in , . (7) 

3 
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Correction to the saddle point approximation is higher order in the small parameter 1/ z n . Examples of diagrams of 
0{\/ z n ) is shown in Fig. 1(b). In addition to the self-energy correction to the fermion level, vertex correction appears 
which contributes to the local field beyond the MFA. Note that to this order the polarization function is diagonal in 
site indices. Only from the order of 1/z^, there emerge off-diagonal contributions. An example is shown in Fig. 1(c). 
Thus we see that in the limit of large dimensions the MFA becomes exact, and that the polarization function H is 
diagonal in site indices to the next leading order. The renormalized exchange interaction J^- is given by 

(^-% = (^%-IW (8) 

Once the renormalized exchange is determined, the (differential) magnetic susceptibility is derived as the two-particle 
Green function of the fictitious fermions. Namely the (renormalized) cumulant average is related to the susceptibility 
as 

Xij = Wlfif}fj)c = P(SaMj), (9) 

where (• • -) c means the cumulant average, and <5<7; = Ui — (<7j). In terms of the polarization function the susceptibility 
satisfies the Dyson-type equation: 

Xij = U l (S lJ + JilXlj)- (10) 
l 

This equation is also written in the matrix form as 

x-^ir 1 -./. (ii) 

In the next subsection we derive the susceptibility explicitly. 



B. optimization of polarization function 

In a theory accurate only up to 0(1/ z n ), there is some ambiguity how to include higher-order terms. We respect 
the self-consistency or, equivalently, the variational property of the theory in choosing higher-order terms. JTT| ] Let us 
consider the perturbation expansion of the thermodynamic potential f2 = —ThxZ. In order to make a renormalized 
expansion, we introduce a fictitious field Uij which shifts (J _1 )y to (J^ 1 ).^ +Uij in H^. With infinitesimal change of 
Uij , the corresponding change of O occurs as 

25n = J2(&h)cSu tJ . (12) 

ij 

Here we consider the case where Suij couples on the average with only the connected part of (<pi<pj), i.e. the cumulant 
average: 

(<pi</>j)c = (<t>i(f>j) - (<t>i)(<f>j) = TJij- 

An example of external fields without coupling to the disconnected part {<pi){4>j) is the case where Suij has a modulation 
corresponding to a finite momentum with hi constant. 

The Dyson equation gives the relationship between and the bare exchange as J -1 = J -1 + u — II in the matrix 
notation. By using 8u = <5J _1 + SU, we can eliminate the fictitious field. Then we get 

2/3SQ = Tn5(ln J" 1 + JTl) - <5${J}, (13) 

where Tr is the trace over the site index and the functional J} is introduced so as to give 5&/5J = II. Equivalently, 
${ J} is written symbolically as 

${J} = Tr£— ^n„{J}J, (14) 

where Tl n is the polarization function made up of all n-th order skeleton diagrams with respect to J. Upon variation 
of J in an n-th order skeleton, we have n ways of choosing J. Hence the denominator is cancelled in the final result, 
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and 5Q/8J in eq.(^4[) indeed gives II. Now we integrate eq.(|l3|) with respect to the exchange interaction ranging from 
J to J. The result is given by 



2f3(Vt - fi ) = Tr [ln(J _1 J) + J ~T l - l] - 3>{ J}, 



(15) 



where 



comes from the lower limit of integration and corresponds to the MFA. We interpret eq.(|l5|) as giving a variational 
form which is minimized by the optimum J. Since II is diagonal in site indices, only the diagonal part of enters 
the expansion eq.(|l^). Furthermore with the optimum choice J -1 = J -1 — IT, we have Tr(JJ _1 — 1) = Tr(j£>IT) 
where Jd is the diagonal part of the matrix J. Thus except for the term ln(J _1 J), only the diagonal part Jd of J is 
relevant to — Qq. 

The foregoing rearrangement guides us how to select terms of higher order for self-consistency of the theory. For 
explicit calculation of O we decompose H^ as 



H{ — Hr> , 



(16) 



where 



Ha 



r, )ij a i 



Hi = 



1 - 



-J i 1 5fa 2 + (mi - cri)S(f)i - (hi + a l )a i 



(17) 

(18) 
(19) 



Here 5 fa = fa — en is the deviation from the average value cti, J t 1 — (Jo) 1 + 11; and m, = ^2j(J 1 )ij a j- The 
quantities Jij and are to be determined variationally. We regard Hq as the unperturbed Hamiltonian. Namely we 
have 



Z = Z G {cx P [-P(J2 Hi - H D )]) G , 



(20) 



where Zq — Tr exp(— (3Hq) and the average is taken with respect to the distribution cxp(— (3Hq). The unperturbed 
thermodynamic potential £}q = — T In Zq is given by 



n G = -Z TtH jj-1) + 1 ^(J-'hm, 



(21) 



where In J -1 enters via the normalization given by eq.(||). By comparing eqs.(|l5|), ( p0| ) and (pi]), we see that the 
average in eq.(20) should rather be taken with use of exp(—f3Hjj) in order to achieve self-consistency up to 0(l/z n ). 
Thus we obtain the partition function as Z — ZqZl/Zo- Here Zl is the product of local contributions given by 



Z L = [](2tt/3J j )- 1 



/2 



J f3dfaexp(-f3H t ) 



(22) 



and Z D = Yli(Ju/ Ji) 1 ^ 2 comes from H D . 

The self-consistent renormalization up to 0(X/z n ) is carried out by requiring that the thermodynamic potential 
of the system be stationary against independent variations of IT and a^. In accordance with the partition function 
we decompose Q as 



n{u,a} = n G + n L - n D , 

where each suffix corresponds to that of the partition function. We write VLl = Qi with 

J7, = —Tin Z;. 



(23) 
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The simplicity of the Ising model allows us to obtain analytically. Namely we first perform the cj> integration and 
then sum over a. The result is given by 

fli = ~Ji(m% + 1) - Tln[2coshj9fti], (24) 
hi = hi +a.i - Jiini. (25) 

We interpret hi as the effective field at site i. Note that Jiirii represents the reaction field |l8[ which is absent in the 
effective field of the MFA. By variation with respect to a,j we obtain 

= ~ J i( J ~%][ m * ~ tanh(/3^)] = 0, (26) 

•* i 

or equivalently 

m l = tanh(/3hi). (27) 

We emphasize that the reaction field enters automatically in our theory as a consequence of the variational principle. 
On the other hand variation of IT gives 

^ = ^[1 - (M- 2 {J%rW^ - (6$)} = 0. (28) 

Here we have used the matrix notation for J 2 and the average (8<t> 2 ) should be taken with reference to Z^. From 
cq.(|2|) we obtain 

(64> 2 )=TJ u =TJ i (l-Il i j i )- 1 . (29) 

Equation (^9|) asserts that IT gives the renormalization of the bare on-site interaction Ji to Jn just as it does the 
renormalization from to . 

We now derive the susceptibility by using the formula: 

on 

Xij + Pmirrij = /3{ai<rj) = -2/3-^—. (30) 

In taking the derivative of il we care only the explicit dependence of because of the stationary property represented 
by eqs.(^) and (^8|). Then the contribution comes only from giving 

Xij = [n(i - Jn)- 1 ^, (31) 

where the matrix notations are used. This result is consistent with the perturbation formula given by eq.(|lO|). 
We quote another relation which is obtained by partial integration of the quantity: 



(32) 



Namely we obtain 

(tr«7j) = ^(J-^iiiMmKJ-^mj - TiJ-%. (33) 

Im 

This formula gives the rigorous relationship between the susceptibility and (<pi4> m )- Substituting {<j)i<j) m ) = TJi m +aia m 
and J -1 = J -1 —II, we obtain the same result for x%j as given by eq. (|3l|) . Thus we see that the single-site optimization 
is consistent with the general property given by eq. ([33]) . 

By using eq.(^) together with eq.(|33|) applied to the effective single-site system, we obtain the local susceptibility 

Xu as 

Xu^^-Ji' (34) 



With matrix notations \ f° r t ne susceptibility, and xl for the local susceptibility, eq.(|31j) is equally written with the 
help of eq. (|34|) as 



G 



(35) 



where J and J are also considered as matrices. In the case of homogeneous magnetic field, the system acquires the 
translational invariance. Then we should recover the local susceptibility by the momentum average of x(l)- The 
self-consistency relation is given by 



1 = Av- 



1 



9 l-(Jq-J) X L 

where Avq — N^ 1 and Jq is the Fourier transform of Jy. One can also write the same relationship as 

J 



(36) 



i - jn 



J q 

Av - y 
q 1 



J rr 



(37) 



Now we mention a drawback of the SMA that the Ward-Takahashi identity (WTI) is violated. In order to see 
a consequence of the violation we derive the susceptibility Xij — drrii/dhj by direct differentiation. We obtain from 
eqs.pl) and © 



Xij = - m?) 



£(■* 



I - o l iJi)xij - mi—— 



(38) 



The local susceptibility is given by Xu — — fnf) independent of the interaction. Hence we recover eq.(|3l|) only if 
we can neglect dJi/dhj. Fortunately this is indeed justified in the zero-field limit since the time reversal invariance 
requires J to be an even function of magnetic field. At finite field, however, the susceptibility depends on whether 
one derives it from fluctuation formula or from the thermodynamic derivative. This drawback of the theory should 
be kept in mind if one discuss the property in finite magnetic fields. Jig ] 

The violation of the WTI originates from the site-diagonal property of the polarization function, and is very hard 
to remedy completely. The WTI requires consistency between quantities with different orders of magnitude in l/z n . 
For example, all the site-diagonal self-energy diagrams in the fermion representation are included in the SMA since 
they are of 0(1) and 0(l/z n ). However, the WTI requires inclusion of polarization diagrams of 0(1/ z^) as well. An 
example is the one shown in Fig. 1(c) which is generated from an 0(l/z n ) part of the self-energy. In the absence of 
magnetization, this diagram and related ones vanish by the particle-hole symmetry. Thus one recovers the WTI in 
this limit. 



C. Entropy and specific heat 



The entropy S of the system is derived by the standard formula S = —dQ/dT. The stationary property of f2 leads 
to enormous simplification; one can neglect the implicit T-dependence of the variational parameters dj and . Hence 
we obtain 



S = -Tr[ln(l - Jn) - ln(l - JU)} + ^ S, 



(39) 



where Si = ln(2 cosh/3/ij) — flhirrii is the contribution of the spin at site i. The form of Si becomes the same as that 
of an isolated spin if one represents Si in terms of m^. Namely we have 



Si = In 2 - ^(1 + m/) ln(l + to.;) - ^(1 - m/) ln(l - m,). 



The specific heat C = T(dS / dT)h of the system is given by 

C 



dJi ( 1 + rrii \ drrii 
Xu^^ + In 



dT 



1-mi) dT 



(40) 



(41) 



where we have used the cancellation implied by eq. 
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In the following we derive the specific heat explicitly in the case without magnetic field. We set m, = assuming 
the paramagnetic phase and remove the site indices of all local quantities. At high temperatures, it is easy to derive 
the leading term as 

J = Av J 2 /T = A/T. 

Then the entropy and specific heat is given by 

S -In 2- A/T 2 , C-A/T 2 . (42) 

On the other hand, for general temperature it is convenient to represent dJ/dT in terms of susceptibilities. The 
temperature dependence of J is derived by the use of the self-consistency relation eq.(|36|). After some algebra we 
obtain a compact formula: 

which is assured to be positive by the inequality Avq(xq) > (Avqxq) 2 - It is evident from above that the specific heat 
and the susceptibility shows the anomaly at the same temperature, and that C in the paramagnetic phase becomes 
larger as intersite correlation develops. This is in sharp contrast to the MFA where C — for T > T c . However, the 
critical property of the present theory is the same as the MFA. 

III. GENERALIZATION TO QUANTUM SPINS 

A. Self-consistent equations for the Heisenberg model 

In this section we extend the formalism developed in §2 to quantum models. As a representative of quantum models 
we consider the generalized Heisenberg model given by 

H=-^J ij S i -S j -J2h-S i , (44) 

ij i 

For clarity we assume in this section that the magnetic field is homogeneous. We emphasize that the formalism 
can be equally applied to systems without the translational invariance and with lower symmetry. The path integral 
representation is accomplished by the coherent-state representation of spins. Equivalently one introduces fictitious 
fermions as 

Si = Y,fl<r<*/3fi0, (45) 

a 13 

where er is the vector composed of the Pauli matrices and the constraint f} a fi a = 1 is imposed. Note that S$ is 
twice the usual spin operator. Then the trace over spin degrees of freedom is accomplished by the path integral |2l[| 



(46) 



where fi a (r) and // q (t) are the Grassmann numbers and is the Lagrange multiplier fields to enforce the constraint. 
Sb stands for the Berry phase term. 

With this preliminary we obtain the expression of the partition function 

VSex V [-S B - d,T H(t)]. (47) 



o 



This form enables us to utilize the Gaussian identity eq.(^|) for each imaginary time interval At, since in expanding 
the exponential we can forget about the non-commuting nature of quantum operators. As generalization from the 
classical case we introduce the time-dependent vector field (/^ (r) for each site, and write Z as 
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Z = J VSV<Pcxp[-S B - drH^T)], (48) 

where (r) is given by 

H f = \ £( J_1 )^ • - + h ) ■ s *- ( 49 ) 

ij i 

Here all time-dependent quantities are specified at time r. 

The renormalization goes parallel to the classical case if we regard Jij as small quantities of 0(1/ z n ). Namely we 
introduce a as the static molecular field at each site. The polarization function now has retardation effect, and Tl(iv n ) 
represents a Fourier component where v n = 2nnT is the Matsubara frequency with n being an integer. As in the 
Ising case, the polarization function is diagonal in site indices up to 0(1/ z n ). However, it has now two components 
U^(iv n ) and H^(iv n ) in the presence of nonzero m. Accordingly the effective exchange interaction Jij(iv n ) is given 
in the momentum space by 

Jq^n)' 1 = Jq 1 ~ n(ii/„), (50) 

which is understood as a 3 x 3 matrix equation with the diagonal matrix n = diag(n ± , n^, n"). Since Jq is a scalar 

in the Heisenberg model, we obtain the diagonal matrix Jq as diag(Jg , Jq , jjj). 

The thermodynamic potential is decomposed as — Qg + — in the same way as in the classical case. Writing 
JD(iv n ) = Ja(wn) wc obtain 

fi G = -|X> In [Jq^/Jg] +I^(J- 1 ) y .a 2 , (51) 

qn ij 

n D = —N^ixHMwnWwn)- 1 ], (52) 

n 

where tr is the trace over the Cartesian indices, and J(w„) _1 = Jo^n) -1 + n(w n ). The nontrivial part is Ql = 
—NTlnZi = Nili where the partition function Z\ of the effective impurity is given by 

Z x = f DSiD^expt- f \tL x (t) - S B \]. (53) 



Here Sbi is the Berry phase term of the site and 

£ i( T ) = 4 f dr'S^r) ■ J(t - r')- 1 ^(r') + [m - S(r)] ■ 5<j>(r) - (h + a) ■ S(r). (54) 
l P Jo 

In the above we have suppressed the obvious spatial index 1, and 

J(t) = ^2 J(iv n ) exp(-ii/ n r). (55) 

n 

One has to solve the single-site problem explicitly to obtain the local susceptibility matrix xl(w„) for given J(iv n ). 
By symmetry these quantities are diagonal with two independent components indexed by _L and ||. In contrast to the 
Ising model, it is not possible to obtain Z\ analytically. Instead, stationary conditions for fl against variations of a 
and H(iv n ) lead to the following relations: 



/ 

Jo 



dr(6<P(r)) = 0, (5<l> a (r)<% (r')) = T [J D (r - r')] a(j . (56) 



Thus we obtain the set of equations which generalize those for the Ising model. For the two components A =_L, || wc 
obtain 

x x L (iv n ) = n A ^„)[i - JVJnVrOr 1 - Av X x (q,iv n )- (57) 

Here the q-dependent dynamical susceptibility X X (q,iv n ) is given by 



9 



X X {q,iv n ) = n A (ii/„)[l - JqT&iiUn)]- 1 = xli^n) {l - [Jq - J X {w)]x X L {w n )} * ■ (58) 
Thus wc obtain the self-consistency equation 

1 = Av{l - [J q - JVrOki^n)} 1 , (59) 

for each Matsubara frequency and each component A. In §5 we generalize these equations to arbitrary localized 
configurations . 

B. Entropy 

In the quantum case the entropy S of the system is conveniently derived by the formula TS = U — Q where U = (H) 
is the internal energy of the system. Another formula S = —dQ/dT is less convenient in the quantum case because 
one also has to take the derivative of Matsubara frequencies. Let us first express f2i by rearranging the perturbation 
terms as in §2 or in the Fermi liquid theory. p2[ In the present case we regard the term 8<p ■ S as the perturbation. In 
order to perform the renormalized expansion we introduce a fictitious matrix field u{t) which increases J(r — t') 
in eq.(||) to J(r - r') _1 + u(t - r'). We obtain 

25ni = J2( 5 ^^Mr'))5u a p(T-r , )=Tj2MT-T%P^(T-T'), (60) 

a/3 a/3 

where Jd(t — t') is generalized from its original meaning with u = 0. The Dyson equation represents the relationship 
in the general case as J^, 1 = J -1 + u — II, and we obtain Su = SJ^ 1 + SH. Then with some manipulation similar to 
the one used in §2, we can integrate eq. (|60|) from the decoupled limit Jd — J to the actual value of Jd- The result is 

2f>!{ J D } = T^tr {- hx[J D {w n )J{iv n )- 1 ] + J D {iv n )J{iv n )~ l - l} - J D } + 2O , (61) 

n 

where 

O = -Tln(2cosh/3|ft, + a\), 

is the MFA contribution. The functional &i{Jd} satisfies the relation 6<f>x/8JD(iVn) = n(w„). As in the classical 
case $i{ Jd} is written symbolically as 

$i{J D } = trY,—-IL n {J D }J D , (62) 

where II„ is a part of the polarization function made up of all n-th order skeleton diagrams with respect to Jd, and 
the frequency summation is implicit. 

We note that m and a are kept fixed during the variation, although actual change of the Hamiltonian does change 
the equilibrium magnetization. This is because the procedure of variation is taken only to use the topological structure 
of perturbation processes. We should also notice the stationary property SQi/SJd = at J^ 1 = J -1 — II in eq. (|6ll) . 

The internal energy is given by 

U = --trJ2 Jq X (q, w n ) - ]-m 2 ^ J tj -Nh m. (63) 



2 

qn ij 



Then from U — f2 we obtain 



S = ~\j2 tT { ln tl - JqTXfiVn)] + 2J D {iv n )Yl{iv n )} + y$i + S , (64) 
q n 

where Sq/N = —f3(h + a) ■ m + ln(2 cosh/3|/i + a\). In contrast to the Ising case, Sq does not have the reaction-field 
correction -Jm. The information of the reaction field is hidden in <&i. Thus without solving the single-site problem 
explicitly, there is not much we can say generally. 
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Nevertheless the result eq.(|54|) is convenient to discuss the limiting behavior. If the magnetic field is absent in the 
paramagnetic phase, we always have So — iVm2. In the high temperature limit, this term alone survives as it should. 
As temperature T decreases, the first term in cq.(^34|) makes positive contribution from extended spin fluctuations. On 
the other hand the following terms up to N<S>±/2 gives a negative contribution to correct the overcounting of on-site 
spin fluctuations. If the entropy given by eq.(p4) does not tend to vanish as T goes to zero, the paramagnetic state 
becomes unstable against a magnetic order. This is always the case for classical spins. The possibility to have the 
paramagnetic ground state in the extreme quantum case is discussed in more detail later. 



C. Elimination of auxiliary fields 



Instead of performing the linked-cluster expansion as described above, it is also possible to integrate away the <j> 
fields in eq.(|54). The resultant Lagrangian C s is given by (apart from the Berry phase term) 

1 rf> i 
£ s{t) = --= / dr'S(r) • J(r - r')S(r') - -m • J(0)m - S(r) ■ (h + a - J(O)m), (65) 

where the matrix j(0) denotes the static component of J{iv n ). Interestingly this Lagrangian contains explicitly the 
reaction field in contrast to the previous expansion. In spite of its simple appearance, treatment of the Lagrangian 
is not easy because of the Berry-phase term. In the next section we show that the equivalent partition function is 
obtained by introducing a fictitious Hamiltonian. 



D. Static approximation 

If the temperature is much higher than T c , only the static component v n = of the Matsubara frequency remains 
important. Then it is reasonable to neglect all the other components. This is called the static approximation. The 
static approximation in the present theory still keeps the non-commuting character of quantum spins. In this respect 
it is different from our previous theory, p4j which replaced the spin vector composed of the Pauli matrices by a 
classical vector. 

In the static approximation, the quantum spin sees a static external field which has the Gaussian distribution 
specified by J. The path integral over S can be replaced by the trace of the density operator with the quantization 
axis taken in the direction of h + a + T£ with £ = j38<p. We then have 

Z x = det(27r/3 J)' 112 J d £ ex P ' J~ l $-m-$j 2cosh/3|/i + a + T£|. (66) 

In contrast to the case of the Ising model, integration over £ does not lead to concise expression. However, the result 
simplifies in some limiting cases. First, by replacing the vector quantities by corresponding scalars, eq.(|66|) is reduced 
to 

Z x = exp[^J(l + m 2 )]2cosh[/?/i], (67) 

where h = h + a — Jm is the effective field with account of the reaction field. Thus the result for the Ising model is 
reproduced. 

On the other hand, if one neglects d/dr in the Berry phase term, the non-commuting nature of spin operators is 
lost. Then we are left with spins behaving as classical vectors. This approximation is justified in the limit of large 
spin If we further neglect the anisotropy in J, we obtain from eq. (p5|) 

Z x = J dfl s exp[^J(S 2 + to 2 ) + PS ■ h], (68) 

where h = h + a — Jm, and the integration is over the solid angle of S. Thus we obtain the partition function of 
classical dipoles with account of the reaction field. The local susceptibility in this case is obtained as 

X [ = PS 2 [x~ 2 - sinlT 2 2;], xt = PS 2 [x^ cothx - x"\ (69) 

where || (_L) is the component parallel (perpendicular) to m and x = [3hS. 
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IV. MAPPING TO SPIN-BOSON HAMILTONIAN 



For numerical calculation, it is often convenient to work with a hypothetical Hamiltonian of the impurity spin 
interacting with bosons representing the effective medium. We shall show that the equivalent Hamiltonian is given by 



H = }Z 



qX L 



u q \b\ x b qX + 



fJ 



^JNuj q x 



(m - S) ■ e qX (b qX + 6j x ) 



-(h + a)-S. 



(70) 



This form is reminiscent of the spin-boson model investigated in another context. pq-|28|] The same structure of 
perturbation diagrams is obtained whichever the boson field or the <p field is used. Let the partition function of the 
impurity plus bosons be given by Zh- Then the partition function Z\ of the effective impurity is given by Z\ = Zh/Zj, 
where Zj, is the partition function of bosons without the impurity. 

In order to characterize the boson system leading to the equivalent partition function we replace 84>{t) by boson 
operators as 



S(f>(r) 



V = 



b g x(r) 



(71) 



where b q \(r) is the annihilation operator of a boson with polarization e q \, and g is the coupling constant. The role 
of g is only to adjust the dimension, and hence can be taken unity in practical calculation. Without loss of generality 
we can take the one-dimensional boson system with q > 0. The correlation function is replaced by the bare boson 
Green function D\(q,ii/ n ). 



N 



^ qX qX D x (q, iv n ) exp(-w„r), 



U qX 



(72) 



where D\{q, iv n ) — (iv n — <jj q \)~ x — (iv n + uj^)" 1 with to q \ being the boson frequency. We write A as _L (two-fold 
degenerate) and || according to the anisotropy introduced by the finite magnetization. Then the bare spectral intensity 
is given by 



ilmJ» = £ 

TT ^ ' 



N0J qX 



[5{u 



J qX) 



8{u- 



->qX 



= 9 2 a-L ( dq 

2-KUJ q \ \dbJqX 



(73) 



with q satisfying uj q \ = lo and ah the lattice constant in the rightmost expression. Here we have assumed that uj q \ 
increases monotonically with q. Thus for given ImJ(w) and g, the spectrum uj q \ is derived from eq.([73]). For example, 
if the boson has a spectrum uj q x oc q i ' p for small q, we obtain ImJ A (w) oc lu p ~ 2 for small u. 

The ground state of the impurity without magnetic field becomes either doublet or singlet depending on the details 
of J(t). In the former case, the self-consistent solution actually drives the ground state to magnetic ordering. On the 
contrary, if the effective impurity has the singlet ground state, the system as a whole also has the singlet ground state. 
This singlet state is reminiscent of the RVB picture. ]3l| Technically eq. ( f70| ) is not exactly the same as the standard 
spin-boson model which is anisotropic in the spin space. In the bosonization approach to the Kondo problem, |?2 33 
one ends up with an equivalent classical two-component gas in one dimension with the interaction decaying as inverse 
square of the distance. In our case, the long-time behavior is given by J(r) oc r -2 if ImJ(w) cx w for small frequencies. 
However, the resultant model is not an Ising model as seen from eq.(pq). It is possible to convert the model to an 
Ising-like one by eliminating the interaction part with S z by a canonical transformation. fl32[ Then dS z /dr plays the 
role of instantons. It remains to see how the self-consistency condition determines the actual shape of J(r). 

Among appropriate methods to solve the effective impurity problem, we mention numerical methods such as the 
numerical renormalization group, Q| the quantum Monte Carlo |l^j3^], or the resolvent method |lj] all of which 
respect the strong on-site correlation. It seems that the most convenient way to obtain the self-consistent solution is 
to use the numerical iteration. The iterative procedure starts from a trial ImJ(w) and the resultant J(iu n ). Then 
explicit solution of the Hamiltonian gives Xi(w„), and hence Tl(iv n ). Substituting these single-site quantities to the 
self-consistency equation (59) gives the second trial J{iv n ) by 



J{iv n ) = n(w„ 



|a v [ii(w„)- 



JqY 



(74) 



Then one continues the second step of the iteration. 
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V. APPLICATION TO HIGHER MULTIPOLES 



The present formalism is flexible enough to be applied to various systems of localized electrons with two-body 
interactions. Let a local electronic configuration a at site i be represented by \i, a). Then the transition to this 
state from another configuration \i,b) is described by the Hubbard operator Xf b = \i,a)(i,b\. In the case of a = b 
the X-operator is equivalent to the projection operator onto the state a. For notational simplicity, we work with 
hermitian operators xj ab} = {Xf + X$ a )/2, x\ ab] = (Xf - X ba )/(2i) and write them as X? where a represents 
either {ab} or [ab]. Then we consider the Hamiltonian given by 

H = -\ E E j ^ x " x j + E e ^ aQ , (75) 

ij a/3 ia 

where the second term represents energy levels with possible splittings caused by magnetic field and/or by crystalline 
electric field (CEF). The exchange interaction has a symmetry much lower than the point-group symmetry at each 
site. Thus it has many nonzero elements with off-diagonal indices in general. 

For the partition function we use the path integral over Xf which can again be performed with use of fictitious 
fermion or boson operators. As in the spin case we impose the constraints that the sum of the projection operators 
be unity at each site. The fluctuating field <f>° has the same set of indices as the X-operators. Then the partition 
function is given by 



Z = / VXV(/)exp 



■S B - drH^r) 
Jo 



(76) 



where Sb is the Berry-phase term and 



\ E E( J-1 )yW + E + E to*™ ( ?? ) 



* 2 

ij a/3 



Here J 1 is the inverse of the double-indexed matrix composed of J"j ■ One can follow the same procedure as in §4 
to renormalize the Hamiltonian up to 0(l/z n ). 

We define the Green function of the fluctuating field by 

Tjf(r - r') = (^(r)^(r'))c = (C(r)^(r')) - <^(t)><<^(t')>, 

where the average is to be taken over the distribution specified self-consistently. To derive this we introduce the 
polarization matrix J\(iv n ) by J -1 = J -1 — n where all quantities are matrices. The important point to note is that 
n Q (r) is diagonal both in site and a. The latter follows from the point group symmetry. By the same reason the 
matrix Jjj is also diagonal with respect to the index a. 

The renormalization is carried out by requiring that matrix of the renormalized single-site propagator in Z\ be 
equal to Jd- The self-consistency equation is given by 

1 = Jj E i 1 - t J 9 - -H^XL^y 1 , (78) 
Q 

where 1 in the left hand side is the unit matrix. The matrix Jq is composed of the Fourier transform of {J^f*}, 
and XL(iv n ) is the local dynamical susceptibility matrix. The effective interaction matrix J[iv n ) satisfies the relation 
X^ 1 = n _1 — J. The full dynamical susceptibility matrix x(<jf, iu n ) is given by 

X(q>„) = n(H/„)[l - JqTliiUn)]- 1 . (79) 

These relations are just generalization of those derived in previous sections. 

The matrix equation is greatly simplified if one uses the basis set of irreducible representations. We sketch the 
procedure assuming the spherical symmetry around the effective impurity. This is merely to make notations familiar, 
and can be straightforwardly generalized to the point-group symmetry. p9| The irreducible tensor operators X(lm) 
are given by 

X(lm) = J2(2J + iy 1/2 (Ja\lmJb)X ab , (80) 
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where (Ja\lmJb) is the Clebsch-Gordan coefficient. Then the local susceptibility matrix in this basis set has only 
diagonal elements as given by 



{lm\xL{iVn)\l'm') = 5w5 mm ,x ( L ( i,y n)- 



(81) 



Since I runs from to 2J, we have (2J + 1) independent components for XL(iv n ). The same number of components 
are present also for U(iv n ) and J(iv n ). Thus eq.(|78]) is reduced to (2J + 1) scalar equations. 

It is convenient to introduce the Hamiltonian of an impurity coupled with fictitious bosons as in the case of the 
Heisenberg model. The Hamiltonian without symmetry breaking fields is given by 



qX 



u q \b\ x b qX + 



y/Nu q x 



(82) 



where the mode A represents the set (Im) of the irreducible tensor. We note that the phase transition to the ordered 
phase is signaled by 



det[l - J q U(0)] = 0, 



(83) 



which reduces to the MFA if one uses the zero-th order result for the polarization function. Namely all the fluctuation 
effect is encoded in IT(0) in the present theory. As temperature decreases, the determinant becomes zero first for the 
most favorable order. If Jq can be regarded as a scalar, there is no coupling among different irreducible representations. 
The opposite extreme case is that q is situated at a low symmetry point of the Brillouin zone. Then many components 
of 11(0) couple in eq. (|S3|) . 



VI. DISCUSSION 



We compare the present theory with previous effective medium theories for fermions. Both theories have similar 
sets of self-consistent equations. Similarity to the fermion theory is seen by the correspondence J G and II £ 
where G is the single-particle Green function and £ the self-energy of fermions. The intersite exchange interaction 
in fermion models is generated by virtual particle-hole pair excitations. Hence a loop of exchange interactions, which 
is taken into account here in deriving the polarization function, is equivalent to a loop of particle-hole bubbles in the 
fermionic theory. However, the intersite processes so far included in the fermionic theory just corresponds to tree 
diagrams of particle-hole bubbles. Thus for magnetic properties the fermionic theory is still in the mean-field level. In 
order to study the itinerant magnetism beyond the mean-field theory, it seems more convenient to use the t-J model 
instead of the Hubbard model. Then the single-site optimization of not only the polarization function but the fermion 
self-energy is accomplished by the use of Grassmann auxiliary fields in addition to the <j> fields. This subject will be 
treated in a separate paper. 

Concerning the possible Kondo-type effect, the Heisenberg model with only the nearest-neighbor interaction J in 
the unfrustrated lattice is not a good candidate. As compared to the characteristic scale z n J for the Neel state, the 
singlet cannot take advantage of many neighbors and the characteristic energy is J. Hence the ground state of the 
nearest-neighbor model becomes the Neel state in high dimensions. On the contrary, if there is substantial frustration 
in the longer-ranged exchange interaction, the critical temperature for the magnetic ordering may be suppressed 
severely. The Kondo temperature is determined by the easiness to form singlet pairs between many neighbors. Hence 
presence of the frustration should work differently. As the simplest example leading to the singlet ground state we 
consider the case where Jy is a negative constant —K independent of the distance. Then the Heisenberg model can 
be solved exactly since the energy is given by 

E = 2KS 2 - C, (84) 

where S = V\ Si/2 is the total spin of the system and C = (3/2)KN. Thus any singlet state gives E — —C as 
the ground state energy. Slight modification of the range of exchange interaction lifts the degeneracy among various 
singlets. It appears possible that certain modification keeps one of the singlet states as the ground state. 

Another important parameter to control the competition between different ground states is the number of compo- 
nents. This number n increases from 1 for the Ising model to 3 for the Heisenberg model. The number 3 is interpreted 
as 2 2 — 1 where 2 is the spin degeneracy and —1 is to remove the identity operator. In the case of four- fold degenerate 
CEF states as in CeBg , we have n = 4 2 - 1 = 15. §6) In the case of TmTe, on the other hand, the CEF splitting seems 
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to be very small as judged from the observed Curie constant [Q which is very close to that of Tm 2+ , and the low- 
energy features observed in the inelastic neutron scattering. [j37| If we assume the degeneracy of the Hund-rule ground 
state with J = 7/2, the number of components is as large as n — 8 2 — 1 = 63. Although the susceptibility matrix 
breaks up into its irreducible components, the overall number of components is still a crucial parameter because of the 
constraint X aa = 1 . In a forthcoming paper we shall report on numerical results of the static approximation for 
the quadrupole order. It is found there that the transition temperature to the quadrupole order is more strongly 
reduced from the MFA value as the number of components increases. 

In summary, we have presented a self-consistent dynamical theory for quantum spins and multipoles. The theory 
is a natural generalization of the spherical model approximation for classical models. Detailed dynamical results 
obtained by numerial calculation will be presented in future papers. 
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